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FORMULAS FOR CONTINUED FRACTIONS 
AN AUTOMATED GUESS AND PROVE APPROACH 

SEBASTIEN MAULAT AND BRUNO SALVY 


Abstract. We describe a simple method that produces automatically closed 
forms for the coefficients of continued fractions expansions of a large num¬ 
ber of special functions. The function is specified by a non-linear differential 
equation and initial conditions. This is used to generate the first few coef¬ 
ficients and from there a conjectured formula. This formula is then proved 
automatically thanks to a linear recurrence satisfied by some remainder terms. 
Extensive experiments show that this simple approach and its straightforward 
generalization to difference and g-difference equations capture a large part of 
the formulas in the literature on continued fractions. 


1. Introduction 


Continued fractions are well known for their approximation properties, their use 
in acceleration of convergence and analytic continuation, as well as their applica¬ 
tion in proofs of irrationality. Any formal power series can be converted into a 
corresponding continued fraction (C-fraction) 
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classically denoted oo -I- fUnW ai(z), 02 ( 2 :), • • • ], where oq is a constant 

and ai{z) are nonconstant monomials for i > 0 that are called partial numerators. 
In the frequent case when all the exponents are equal to 1, the C-fraction is called 
regular. Truncating a continued fraction after its nth term gives a rational function 
which is called its nth convergent. There is a one-to-one correspondence between 
power series and C-fractions. It is easily computed by a sequence of extractions of 
the constant coefficient, division by the variable and inversion. This conversion is 
available in the major computer algebra systems. 

In several isolated cases, the coefficients am are known to possess a closed form, 
as in the following formula for exp(z): 
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or more compactly 

(2) ai=z, a 2 k =-z/{2{2k-1)), a 2 k+i = z/(2{2k + 1)). 

Such formulas are the object of this work. A number of them are listed in the 
classical handbook by Abramowitz and Stegun [1], or in its successor [15] and the 
most extensive list to date is the recent handbook by Cuyt et alii [10]. Our aim is 
to derive many of these formulas automatically, starting from a description of the 
function to be expanded in continued fraction. 

We concentrate on functions that are given as solutions of ordinary differential 
equations with initial conditions (or difference or g-difference equations, see §6). 
Our approach can be summarized as follows. First, the differential equation and 
initial conditions are used to generate the first terms of the power series expansion 
of the function. This power series is then converted into a continued fraction. The 
coefficients of this continued fraction are then “guessed” by variants of rational 
interpolation. When this guessing phase is successful, a new power series is defined 
by this guessed continued fraction expansion. It remains to show that this power 
series satisfies the differential equation (the initial conditions being correct by con¬ 
struction). The key point in this proof is Theorem 4, stating that the (properly 
normalized) evaluations of the differential equation on the successive convergents 
to the continued fraction satisfy a linear recurrence, that can be computed. In 
all cases, after an operation we call “reduction of order”, this recurrence exhibits 
a growth in the valuations that is sufficient to conclude the proof. A surprisingly 
large proportion of known explicit continued fractions are thus obtained completely 
automatically. 

Classically, a very effective method due to Gauss derives formulas for continued 
fractions starting from the contiguity of hypergeometric series. Specialization of 
the parameters then leads to formulas for elementary or special functions [12, §6.1]. 
This leads to explicit continued fraction expansions by recognizing the function to 
be expanded as a special case of a quotient of contiguous hypergeometric functions 
and then relying on a small table of such explicit formulas. These quotients satisfy 
Riccati equations, so that they are covered by our approach, which is not limited to 
them (see §6) and proves more suited to the targetted application to the Dynamic 
Dictionary of Mathematical Functions [5] . This is an online encyclopedia of special 
functions, where the formulas are all generated by computer algebra algorithms 
from differential equations, in many cases along with a human-readable proof. In 
this context, it makes sense to avoid any table lookup and generate formulas and 
proofs for continued fractions directly from the differential equation. 

The work closer to ours is the investigation by Chudnovsky and Chudnovsky [8] . 
They used computer algebra in the study of formulas for continued fractions. Their 
aim was to classify all functions possessing continued fractions with explicit formulas 
of various types and relating them to Painleve transcendents. In contrast, we focus 
on one function that is given as input, and heuristically produce a rational continued 
fraction expansion when possible. 

This article is structured as follows. Section 2 gives an overview of our method 
on the example of the tangent function. Next, Section 3 presents a heuristic of 
independent interest that reduces the order of a recurrence given initial conditions. 
This plays a crucial role in the proving phase of our method. Section 4 is a brief 
account of what guessing means in this context, while Section 5 is the heart of this 
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work and shows how proofs are achieved automatically. Finally, Section 6 presents 
experiments with this approach. 

2. Detailed Example: tan 

The tangent function can be defined by the Riccati equation 

(3) y' = i + y\ y(0) = 0. 

The first 15 coefficients of the unique power series solution are easily computed from 
the differential equation (see Proposition 3 below for existence and uniqueness). A 
conversion into a continued fraction gives the coefhcients 

[0,2, -zV3, -zVl43]. 

The general formula can be deduced from these first terms by rational interpolation, 
which leads automatically to the (so far conjectural) formula 

(4) ai{z) := z; a„( 2 ) := —z^/((2n — 3)(2n — 1)), n > 1. 

Next, we turn to the automatic proof of this formula. The strategy is to prove 
that the sequence of rational functions defined by truncating (4) after the nth term 
for n = 1,2,... converges to the formal power series solution to the differential 
equation (3). More precisely, let /„ be defined by 

/„ = ^ := [0, z, -zV3, ■ ■ ■, -zV((2n - 3)(2n - 1))], 

where the rightmost term denotes the finite continued fraction. Then the proof will 
be completed by showing that 

val(tan —fn) ^ oo as n ^ oo, where val denotes the valuation of a power series: 
val (^ CiZ*) := min{z > 0 | Ci ^ 0}, 

i>0 

with the convention val(O) = oo. Proposition 3 below shows that it is sufficient to 
prove that val(T>(/„)) —>■ oo, where 

It is classical that the numerator and denominator of the convergents of a con¬ 
tinued fraction are related to the coefficients a„ through a linear recurrence: 

(P_l,Po) = (1,0), Pn = Pn-l + anPn-2, n > 1, 

(Q — I 7 — (0, 1), Qn — Qn—1 “f ^nQn—2i U ^ 1. 

In view of (4), it follows that for all n > 0, <5^(0) = 1. Thus, the valuation 
of I){fn) is that of its numerator 

(6) := -Ql- P^ - PuQ'n- 

Using (5) to rewrite Pn+k and Qn+k in terms of Pn, Pn+i, Qn, Qn+i, it follows that 
any shift Pln+k [k G {0,1,2,...}) can be rewritten as a linear combination of 

Pn+iQn+j , Pn+iQn+jt Pn+iPn+j , Qn+iQn+j , 

for i and j in {0,1}. There are finitely many such terms, which implies that a 
linear dependency between i/„, ... (ie, a linear recurrence for iJ„) can be 
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computed directly from (5) by linear algebra. This computation produces a linear 
recurrence of order 4: 


(7) (2n + 7)z®i7„ - z\2n + 7)(2n + 3f{2n + 

+ 2z^{2n + 5)(2n + ?,f{2n + + 20n + 21)i7„+2 

- {2n + 5f(2n + lf(2n + 7)‘^{2n + 3fH„+3 

+ {2n + 5f{2n + lf{2n + 7f{2n + 2,fHn+i = 0. 

This recurrence is satisfied by all sequences defined by (6), with and Q„ arbitrary 
solutions of (5). Using the actual sequences and Qn provided by the continued 
fraction gives the first values of 


’ ’ 9’ 225’ 11025’ 893025' 

From there, automatic guessing again suggests the following simpler recurrence 
for iJ„: 

(8) (2n + l)2i7„+i - = 0. 

And again, this recurrence admits of an automatic proof: the right Euclidean divi¬ 
sion of the fourth order recurrence operator from (7) by this first order one has a 
remainder equal to 0. This shows that the solution of (8) with the initial conditions 
given above coincides with the solution of (7) with the same initial conditions, and 
thus the numerator of P(/n) satisfies (8). On this last recurrence, the increase of 
the valuation with n is clear and this concludes the proof that /„ converges to tan 
and thus that the power series defined by the continued fraction (4) is that of tan. 

In summary, starting with the differential equation (3), this method produces 
and proves automatically the general term of the famous continued fraction 


tan z =- 

1 - 


zV3 


z^/((2n — 3)(2n — 1)) 


that was the basis of Lambert’s proof that tt is irrational. 

3. Reduction of order by Guess and Prove 

The transformation of the large recurrence (7) into the shorter one (8) makes 
it possible to prove automatically that the valuations vali7„(z) increase with n. 
This transformation turns out to play a role in most of the examples dealt with 
in our experiments. It is actually of more general interest: the closure properties 
enjoyed by the class of D-finite series or P-recursive sequences give rise to operators 
satisfied by products or sums of zeroes of such operators [17, 16]. These operators 
annihilate all possible cases and are potentially of large size, while operators of 
smaller order may exist for the specific solution of interest. Such an operator may 
be a right factor of the large one and could be searched for by factoring, but this 
is made difficult by the potentially infinite number of distinct factorizations [18]. 
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Sequences. Let A he a, recurrence operator with polynomial coefficients in n, 
of order denoted by ord^ and leading coefficient lc(^)(n). A sequence {un)n>o, 
abbreviated (un), is said to be defined by the operator A and the initial condi¬ 
tions K. = (ui)jgj, when the value Un is given by /C for n G I and by the recurrence 
operator evaluated at n — ord A otherwise. Note that the set I must contain 

{0,..., ord A — 1} U {z G N I lc(A)(z — ord A) = 0}. 

Algorithm. We now detail an efficient heuristic approach finding such right 
factors, whose complexity is controlled with respect to the size of the large operator. 
The idea is to exploit the initial values of the sequence by a “guess and prove” 
approach. This is described in Algorithm 1. This algorithm takes as input a linear 
recurrence operator A and initial values, as well as an upper bound N on the 
number of coefficients used to find a right factor. It is described here in the case of 
recurrence operators; similar variants apply to differential or q-difference cases. 

The search for a smaller order operator is performed in two main steps, “guess¬ 
ing” and “proving”. First, the input recurrence of order M and its initial values are 
used to compute the first N terms of the sequence. Next, these N terms are used 
to “guess” a linear recurrence. This is done by linear algebra: we search succes¬ 
sively for the existence of a linear recurrence operator Q of order 1, 2,..., M with 
polynomial coefficients of degrees such that the sum of the numbers of undeter¬ 
mined coefficients of the recurrence is smaller than N. The structure of this linear 
algebra problem is exploited by computing matrix rational interpolants [4] (in the 
differential case, Hermite-Pade approximants are used [3]). 

When N is sufficiently large, this linear algebra phase is always successful, since 
it can reconstruct A. The next step is to prove that the recurrence Q obtained from 
the first N terms of (u„) defines the same sequence for all n. The operator Q is 
not necessarily a right factor of A, but could be merely a left multiple of such a 
right factor, the factor itself being too large to be found with N terms only. This 
is related to the typical shape of the order-degree curve [7]. Thus the algorithm 
next computes the greatest common right divisor of Q and A and its numerator 7?., 
obtained by left-multiplication with the least common multiple of the denominators 
of the coefficients. 

At this stage, the algorithm has produced a right factor TZ of A. It is then 
associated initial conditions with which TZ defines a sequence . We 

now prove that if for n G -I- {1,..., ordA — ordT?.}, then Vn = Un for all 

n. The induction on n is based on the following. 

Lemma 1. Ifun = Vn forn < z-|-ord A—I and i+j ^ J for all j G {ordT^,..., ord A} 
then ZXj-i-ordA — T^-i-ordA- 

Proof. The sequences {S^TZ ■ Vn}ordTZ<j<ordA all cancel at n = i. The application 
of A at is a linear combination of them with coefficients that are finite, as 

shown in the next lemma, so that A ■ is 0 at n = z. □ 

Lemma 2. Let A and C be recurrence operators with polynomial coefficients, sat¬ 
isfying A = BC where B has rational coefficients. Then the denominator denB 
satisfies: 

den{B)-\0) C lc(C)-i(0) -f {0, -I,..., -ordS} 
where addition denotes the sumset. 
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Algorithm 1 Reduction of Order 

Input: (A, (wn)„gi) defining , and A > 0. 

Output: {TZ, defining (Mn)„>o s.t. ordTZ < ord A. 

n ^ {un)n=o,...,N-i, Computed using A and 
0 -ir- guessrec(A) 
if ^ ^ FAIL then 

TZ numer (gcdrightC^, G)) 

J (lc(7^)-H0) -ord7^) ON 

^ (^ri)n£ J'+{l,...,ord A—ord TJ} > USing {TZ,{Un)^^j) 

w ^ (Mn)riej7+{l....,ordA-ord-R}, USing {A, {Un)„^x) if = V 

then return {TZ, {un)n^j) end if 
end if 

return {A,{un)^^x) 


Proof. This is seen by following the steps of a right Euclidean division. □ 

In practice, this algorithm is run for increasing values of N = 4,8,16,... and 
stopped when either a factor is found or N is larger than the number of coefficients 
of A. Note however that if A is not irreducible, then increasing N further is bound 
to find a nontrivial right factor. 

4. Guessing Continued Fractions 

The first step of our approach to continued fractions is the automatic discovery 
of formulas for the partial numerators Ofc. This section is very short since this part 
of the computation is straightforward. 

Starting from the differential equation, a first method would be to produce the 
first terms of the series expansion of the function, convert them into the first terms 
of the continued fraction and then use the method of the previous section to look 
for a linear recurrence of size bounded by the number of terms that have been com¬ 
puted. It turns out that in most of the known examples, explicit formulas are of 
rational form (see Section 6). We therefore concentrate on rational coefficients, or 
on interlacing of rational coefficients as in the example of the exponential function. 
This means that the “guessing” stage of our approach relies simply on rational in¬ 
terpolation, problem for which efficient algorithms are known through its relation 
to the extended Euclidean algorithm [11, §5.7]. Moreover, the degrees of the nu¬ 
merator and denominator are generally low, so that a few terms of the expansion 
are sufficient for the computation. 

5. Proving Continued Eractions 

FOR SOLUTIONS OF ORDINARY 
DIFFERENTIAL EQUATIONS 

The proving phase is the heart of our work. It is also where most of the compu¬ 
tational work takes place. We consider first-order non-linear differential equations 
with rational coefficients, ie, y' = p{y), with p G Q( 2 :)[F] of degree d. In particular, 
the case d < 2 corresponds to the classical Riccati equations that are ubiquitous 
in the study of continued fractions, due to their stability under linear fractional 
transformations of the unknown function [2, 10.7]. Explicit solutions for restricted 
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classes of equations have been provided by Euler and Lagrange and more recently 
by Khovanskii [13]. 

Our procedure goes in the reverse direction. The continued fraction with explicit 
rational coefficients that was found in the previous stage defines a power series. The 
aim is to show that it is a solution of the differential equation. 

5.1. Valuations. The following proposition reduces the proof to that of the ulti¬ 
mate increase of a sequence of integers. 

Proposition 3. Let F € K[[V, V]] he a formal power series with coefficients in a 
field K and let {fn{X)) be a sequence of power series in IK[[V]]. Then the differential 
equation Y' = E(V, Y) with initial condition V(0) = 0 admits a unique power series 
solution S{X). Moreover, the sequence {fn{X)) converges to S(X) (ieval{fn — S) 
oo) if and only *//ra(0) = 0 for n sufficiently large andval{ff{X)—F{X, fn{X))) —>• 
oo. 


Note that equations with an initial condition V(0) = a ^ 0 can often be brought 
to this setting by changing the unknown function into a + Y. 


Proof. Recall that the algebra of power series is a metric space for the distance 
induced by the valuation: if / and g are two power series, then d{f, g) = 
where val denotes the valuation (this distance does not derive from a norm). It is 
a simple consequence of the dehnition that Cauchy sequences for this distance 
converge in IK[[V]]. 

The first part of the proposition is a variant of Cauchy’s theorem, whose proof 
is straightforward thanks to Taylor expansions. In detail, the solutions of Y' = 
F{X,Y) with initial condition y(0) = 0 are the hxed points of the operator Q : 
y i-A / F{X,Y); this operator is a contraction: 


va\{g{Yi) - g{Y 2 )) = val (^J F{X,Yi) - F{X,Y 2 )^ 

= val (^J ||(X, Y 2 ){Yi - Y 2 ) + 0((ri - Fa)")) > vaKVi - ^ 2 ); 

this shows both the existence of a solution (start from F = 0, iterate g and use 
completeness) and its uniqueness. 

Next, if val(/(j - F{X, /„)) = K, while /„(0) = 0 = S'(O), then 

S-fn = ^(0) - /„(0) + J {F{X, S) - F{X, /„)) + 0{x^+^) 

= {gis)-giu)) + o{x^+^). 

The previous inequality with Yi = S and Y 2 = fn shows that the valuation of 
the first term on the right-hand side is larger than that of the left-hand size and 
thus val (S' — fn) > K + 1, which shows that fn ^ S. The converse implication 
follows from the continuity of the map Y 1 -^ Y' — F{X, Y). □ 

This proposition extends to more general equations of the type P{z, y,y',..., y*-"^) = 
0, with natural assumptions on the initial and separant of the equation. 
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5.2. P-recursivity and Convergents. Recall that a sequence is called P-recursive 
when it satisfies a linear recurrence with coefficients that are polynomial in the 
index. P-recursive sequences are closed under sum and product and algorithms 
computing the corresponding recurrences are known and implemented [17, 16]. 

The key to our approach is the following. 

Theorem 4. Let {Pk{z)) and{Qk{z)) be P-recursive sequences of rational functions 
in z and let F € ]K(z)[T] be a polynomial of degree d > 0 in Y. Then the sequence 



satisfies a linear recurrence with coefficients in K[ 2 :, fc]. 

This theorem is used when {Pk) and (Qk) are the sequences of numerators 
and denominators of the continued fraction supposed to converge to a solution 
of y' = F{z,y). Its proof constructs a recurrence for Hk from which the increase 
of valuation will be obtained using the reduction of order of Section 3. This is 
sufficient thanks to Proposition 3 and the observation that val Qk = 0, which will 
follow by induction from Eq. (5) and the fact that val Cfc > 0 in applications to 
C-fractions. 

Again, similar results can be stated for higher order differential equations, but 
they proved unnecessary for the continued fractions dealt with in our experiments. 

Proof. Let M be the order of the recurrence satisfied by {Pk)- Using this recurrence, 
all Pk+i, i G N can be rewritten as linear combinations of Pk+j for j = 0,..., M— 1, 
with coefficients in Q(k,z), while the polynomials rewrite as linear combina¬ 
tions of those same polynomials complemented by Pk+j for j = 0,...,M — 1. The 
same argument applies to the sequence {Qk) and we denote by M' the order of the 
recurrence it satisfies. 

The choice of the exponent of Qk makes Hk a polynomial of degree d in Pk, Qk, 
Pi, and Qfc- Thus all the Hk+i for i £ N can be rewritten as linear combinations of 
monomials of degree d in Pk+i,Qk+j, i = 0, ■ ■ ■, M — 1, j = 0,, M' — 1 and their 
derivatives. These monomials are in finite number N, whence a linear dependency 
between Hk,. ■., Hk+N (ie, a linear recurrence of order at most N satisfied by {Hk)). 
It can be found by linear algebra. □ 

As exemplified by the computation in the example of tan in Section 2, the bound 
on the order on specific examples may not be as large as suggested by this proof. 
Our implementation thus proceeds by increasing the order one by one and looking 
for a linear dependency until one is found. The algorithm is outlined in Algorithm 2. 
Its termination is granted by the theorem. 

We state a simple generalization of this result that could be useful in applications 
to continued fractions: if the partial numerators in the continued fraction expan¬ 
sion (1) are of the form r{k)z^^^'^ with r rational and p polynomial, then again, the 
polynomials Hk defined in the theorem satisfy a linear recurrence, this time with 
coefficients that are polynomials in k, z and a finite number of z^ , with i < degp. 
The proof follows along the same lines. 

5.3. Riccati Equations. The case when the polynomial F of Theorem 4 has de¬ 
gree 2 gives rise to Riccati equations that are ubiquitous in the theory of continued 
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Algorithm 2 Recurrence for {Hk)f,yQ 

Input: linear recurrences Cp and Cq of order bounded by M for (Pfe)j,>g 
and 

Output: a linear recurrence Cp for (^I/c)fe>o 
To(fc) ^Hk 

for i = 1, 2,3,... do 

Ti{k) ^ Ti_i{k + 1) with P/c+m, P'k+M^ Qk+M, Q'k+M rewritten in terms 
of values of these sequences with smaller indices, using Cp, Cq and their 
derivatives. 

if the linear equation + Pi{k) in the unknowns cg,..., Ci-i has 

a solution then 

return Hk+i + Ci_iHk+i-i H-h CgPfc = 0 

end if 
end for 


fractions [2, 10.7]. In this case, the computation of a recurrence of the form pre¬ 
dicted by the theorem can be made explicit in full generality. 

Proposition 5. LetK.’^i ^ solution of the Riccati differential equation Y' = 

F{z,Y) where F is a polynomial in K(z)[y] of degree 2 in Y, let (Pk) and {Qk) be 
sequences obeying the linear recurrence Uk +2 = Uk+i + ak+ 2 {z)'^k, with a'f.{z) 0. 

Let finally Hk be defined by 

Hk = Ql [{PklQk)' - F{z, PklQk )). 

Then the sequence {Hk) satisfies the following linear recurrence of order J): 


-Hk+i + — - 


"-fc+i 


a-k 


Ofe+l + 1 

<+i 

ak{ak 


Hk 

1 ) o,k+i{ak+i + 1 ) 


a 




1 


fc+1 


Hk-i 


a' 


k+l 


^kPk-2 




Hk-3 = 0 . 


The shift of the indices (from Hk+i to Hks) is only for readability. A nice 
property of this recurrence is that its coefficients do not depend on the differential 
equation, but only on the sequence Ofe. This persists for higher degrees: a differential 
equation with a cubic right-hand side leads to a recurrence of order 6 that does not 
depend on the equation. 


Proof. The formula is obtained automatically by the method from the proof of 
Theorem 4, on a differential equation with symbolic coefficients. It could also be 
derived by hand. However, once it is given, it is a simple matter to produce a proof: 
inject the definition of Hk into the recurrence, rewrite all the PCs and QCs using 
the recurrence they satisfy in terms of Pk- 3 , Qk- 3 , Pk- 2 , Qk -2 and collect terms to 
observe that the left-hand side becomes 0. □ 


As an example, setting ak{z) = —z^l{{flk — l)(2fc —3)) recovers Eq. (7) obtained 
for tan. 
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Corollary 6. If the sequences [Pk) and (Qk) satisfy a linear recurrence of the 
form Uk +2 = bk+ 2 iz)uk+i + ak+ 2 (z)uk, then the sequence [Hk) satisfies a fourth- 
order linear recurrence obtained by evaluating that of Prop. 5, replacing oi by ai/bi 
and Ok by ak/(bkbk-i) for k>2. 

Proof. This is a classical transformation of continued fractions. Setting Pk = 
Pk/{bi ■ ■ - bk-ibk) and similarly for Qk and injecting into the recurrence equation 
shows that both sequences {Pk) and [Qk) satisfy 

, Ofc+2 

l^fc+2 — “f T 7 life- 

0fe+20fc+l 

Since Pk/Qk = Pk/Qk, the proposition applies. □ 

5.4. Nonregularity and Periodicities. As the example of the continued fraction 
for the exponential function in Eq. (2) shows, not all common closed forms for 
continued fractions are given by one rational function. However, most C-fractions 
formulas in the literature appear to be “periodic”, in the sense that there exists a 
period ^ > 0, and £ sequences (a^),..., that alternately define the partial 

numerators Ok- ak = a!'^ The case £ = 2 encountered for exp is the most 

common, but higher values also happen (e.g., £ = 4 for ip" , where ip = T'/T is the 
logarithmic derivative of the Gamma function). 

This is not a restriction in our approach, by the following. 

Lemma 7. Given a period £ > 0, a sequence {uk)k>o P-fecursive if and only if 
all its subsequences {uek-kj)k>o P-recursive, for j = 0 .. .£ — 1. 

Proof. This lemma is classical. We give a constructive proof for completeness. 
If the sections {uik+j)k>o P-recursive, then their generating series Sj{z) = 
^^>0 are D-finite, then so is so{z^) -I- zsi{z^) • -I- z^~^Si-i{z^) and 

therefore its sequence of coefficients {uk)k>o P-recursive. Conversely, if {uk)k>o 
is P-recursive, then its generating series 5 ( 2 ;) is D-finite and so is its Hadamard 
product with z^ /{I — z^), then also its quotient by z^ evaluated at z"^P and this is 
precisely the generating series of {ui,k+j)k>c, ■ 

In cases like the exponential function, this lemma implies that the sequence 
of partial numerators {ak)k>o itself satisfies a linear recurrence. With the recur¬ 
rences (5), this alone does not imply that (Pk) and {Qk) are also P-recursive for 
in general, no such closure property exists. For instance, the sequence defined 
by Un := nL fc! satisfies a linear recurrence of order 1 with a coefficient {kl) that 
is P-recursive, but («„) itself is not P-recursive, as can be seen from its asymp¬ 
totic growth that is too fast. The crucial property in our application is that the 
sequences (o^) are rational in k. This allows for the following. 

Lemma 8. Let a^,..., a^fP^ be rational functions of k and z, ak be defined for k > 0 
by Ofc := a!'^ and let the sequences {Pk) and {Qk) be defined by the recur¬ 

rence (5). Then {Pk) and {Qk) are P-recursive sequences. 

Proof. The proof is constructive. By the recurrences (5) and the definition of Ok, 
all Pik+j for j = 1,..., 2£ rewrite as linear combinations of Pik+i and Pik with coef¬ 
ficients that are rational functions of z and k. Thus Pek +j, P^(fc+i)+j and Pi(^k+ 2 )+j 
are linearly dependent for j = 0,...,d—1. 



FORMULAS FOR CONTINUED FRACTIONSAN AUTOMATED GUESS AND PROVE APPROACH 


The same reasoning applies to {Pik+i)k>o for any j £ {0, 1} and shows 

that it is a P-recursive sequence and thefore that so is {Pk) by the previous lemma. 
By construction, (Qfc) satisfies the same recurrence as (Pk)- D 

Example. The special case £ = 2 is important in applications. Starting from 
P 2 k = P 2 k-i + o, 2 kP 2 k -2 and its first two shifts, the linear combination P 2 k +2 + 
P 2 k+i — a, 2 k+iP 2 k gets rid of the terms with odd index, leaving: 

(9) p2fc+2 = (1 + a2fc+l + a2k+2)P2k — CL2ka2k+lP2k-2- 

A similar computation would give a recurrence between the terms with odd index. 

The proof of Lemma 8 leads to an algorithm in two steps: compute a recurrence 
for (Pk) and (Qk) and then appeal to Lemma 7. A simpler and faster computation 
proceeds directly from a recurrence for (Ptk), thanks to the following. 

Proposition 9. Let F(X,Y) be a rational function, that is regular at X = Y = 
0. Let a\, i = !,...,£ be rational functions in X and k with positive valuation 
in X. Let (ak), (Pk) and (Qk) be defined as in the previous lemma and (Hk) as in 
Theorem 4- IfvalHke oo as k ^ oo, then the continued fraction IC^^i ^ is the 
solution ofY' = F(X,Y) with T(0) = 0. 

Here again, other values for Y (0) are obtained by a change of unknown function. 

Proof. Since the denominator of F does not vanish at 0, F admits an expansion 
in power series and thus by Proposition 3, the differential equation Y' = F(X,Y) 
with T(0) =0 possesses a formal power series solution S(X). 

The condition on the valuations of the sequences (a^) makes the continued frac¬ 
tion well-defined, in the sense that the sequence of power series (Pk/Qk) converges 
to a power series G(X). Thus if its subsequence fk = Pkt/Qki converges to S(X), 
then we have G(X) = S(X). Induction from Eq. (5) shows that Qfc(O) = 1 for 
all fc > 0 and that Pfe(O) = 0 (by the positive valuation of a}). This gives /fe(0) = 0 
and valiLfcf = val(/(, — F(X,fk)). Thus by Proposition 3, the sequence (fk) con¬ 
verges to 5'. □ 

Example. The proof for the continued fraction for exp from the introduction goes 
as follows. Starting from the recurrences for (Pk) when k is even and when k is 
odd and proceeding as for Eq. (9) yields 

P2k+2 — P2k + _ 2^^ ^2fc-2) 

which is also satisfied by Q 2 k since this computation does not depend on the initial 
conditions. 

Next, turn to the numerator of the evaluation oi y' — y — 1 at y = P 2 k/Q 2 k, 
namely 

H2k = P^kQ'ik - P2kQ2k - P2kQ2k “ Qlk- 

Using Proposition 5, or directly as in Section 2, leads to a recurrence of order 4, 
on which reduction of order yields H 2 k +2 = —z^L72fc/(4(2/c -|- 1)^), which concludes 
the proof. 
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Algorithm 3 Discovery and proof of continued fractions 

Input: Y' = F(z,Y) with F G K(z)[y]; 

a bound JV on the number of coefficients to guess from; 
a bound L on the period to be found. 

Output: In case of success, an explicit expression for the continued fraction ex¬ 
pansion of the solution such that D(0) = 0. 

Compute the first coefficients oi,..., ojv of the continued fraction expansion of 
the power series solution of Y' = F(z,Y) with y(0) = 0. 
for £ = 1, 2,3,..., L do 

Use rational interpolation to compute aj interpolating the subse¬ 
quences (a£i+j)t, j = 0, ...,£- 1. 
if this has been successful theu 

compute a recurrence TZ for Hki, with Flk defined in Theorem 4 
compute a new recurrence TV from 7^ and the initial conditions for {Hk) 
using Algorithm 1 

if TZ' exhibits the increase of {va\H(k) theu 
returu the rational functions aj 

eud if 
eud if 
eud for 
returu FAIL 


6. Experiments 

An overview of the whole approach is given in Figure 3. In practice, N = 
20 and L = 2 have proved sufficient in our experiments except for one case of 
period 4. For the computation of the first terms of the continued fraction, one can 
either compute a power series expansion first, e.g., by Newton iteration [6], or use 
techniques for continued fraction expansions of solutions of Riccati equations [9] . 

Our main experimental result is the following. 

Empirical Observatiou. All the 53 explicit C-fractions formulas of the com¬ 
pendium by Cuyt et alii [10] can be guessed and proved by our approach and its 
variants below. Among them the vast majority (44) are solutions of Riccati equa¬ 
tions, 2 satisfy q-Riccati equations and the remaining 7 satisfy difference equations. 

We now give more detail on the calculations in the differential case and then 
outline the variants of our method in the g-difference and difference cases. An 
implementation in the differential case is available under the form of a submodule 
gfun:-ContFrac of the package gfun (for versions > 3.70). It can be downloaded 
from our web pages. All the examples of solutions of Riccati equations from [10] 
are provided through the associated help pages. 

6.1. C-fractious from Differeutial Equations. In our experiment, the Riccati 
equations were themselves found by a guessing approach on power series expansions 
to small order (less than 30). Depending on how one decides to define the power 
series from the computer algebra point of view, these Riccati equations can also be 
automatically proved to hold. 
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Gauss’s continued fraction. The classical hypergeometric series is 


2 Fi{a,b]c;z) := ^ 

n>0 


ia)nib)n 
(c)„ n\ ’ 


where (a)n is the Pochhammer symbol (a)„ = a(a+l) • ■ • (o + n —1). Gauss proved 
the following identity 

2 Fi(a,b;c;z) _ ^ 

2 ^ 1 ( 0 , &+1; c + 1; z) m=i 1 

_ (k + b)(k + c — a) (k + a)(k + c — b) 

(2fc + c)(2A: — 1 + c) ’ (2A: + c)(2fc + 1 + c) ’ 

for the quotient of two contiguous hypergeometric series. This is the source of 
many continued fractions for special functions by specialization of the parameters. 
If y = 1 + F is the function on the left-hand side, then elementary properties 

of the 2 F 1 that can be derived from the first order recurrences satisfied by its 

coefficients show that 


cz(z — l)y' = a(c — b)z + {c{a — b)z -I- (?)y + (?y^■ 


This is our starting point. From there, it is easy to compute the first 20 coefficients 
and conjecture the formulas for a 2 k and 02^+1 by rational interpolation. As in 
Eq. (9), a recurrence for even indices follows. From Corollary 6, a linear recurrence 
of order 4 follows for the remainder i? 2 fc) that can be either obtained by hand from 
Proposition 5, or by a generic code that searches for linear dependency. Next, 
reduction of order gives a two-term linear recurrence within a couple of seconds: 


H. 


2 {n + a) in — a + c) {n + b) {n — b + c) 

— 2^ ^ n T -^2(n—3) 


2(n-2) 


(2n -I- cf' {2n + c— 1)^ 


and this concludes the automatic proof. 

More parameters. Khovanskii [13, p. 85] gives an explicit continued fraction 
with 5 parameters for the power series solution of the differential equation 


(1 -I- az)zy' + {15 + iz)y + Sy'^ = ez, y{0) = 0 

(an extra parameter k is obtained by changing z into z^ and adjusting the coefficient 
of y'; we have relabeled the parameters). This contains the equation for Gauss’s 
continued fraction above as a special case. 

From there again rational formulas for a 2 k and 02^+1 are obtained by guessing 
on the first 20 values; a recurrence of order 4 can be found for the remainders Fl 2 k] 
Algorithm 1 reduces it to the conclusive recurrence: 


(2n -I- /3)^(2n + /5 — l)^Ff2(n-2) — 

- {an^ + (a/3 -I- qjn + /5j + 5e){an^ + (a/3 - 7)71 -|- 5e)z‘^H2(n-s)- 


Other examples. We also applied our method to a few functions not men¬ 
tioned by Cuyt et alii [10] and in particular found (and proved) experimentally the 
following nice C-fraction for the Airy function: 


z 




(lA^) 



1+ 1C 

m—2 



02 fc = (6fc - l)z^/8. 


a 2 fc-i-i = (6fc -I- l)z^/8. 
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It also follows from known C-fractions for the divergent 2 -^ 0 . 


6.2. q-analogues. The method used in this article also applies to g-analogues. We 
outline the very simple example of the (/-exponential: 

(1-9)” 




m>0 


(1 - q){l - g2) ... (1 _ gm) 


which satisfies the (/-differential equation 

nn\ eq{qz) - eq{z) 

(10) _ e,(.) = 0. 

The classical exponential is obtained as the limit when (/—>!. The first coefficients 
of the continued fraction expansion let one guess ai = z, 

"'=-1(1-(/)z q^'^{l-q)z 


a2k = —- 


(l + g'=)(l-g2fc-i)’ “2fe+i 

a clear generalization of the continued fraction (2) for exp. In order to prove this 
continued fraction, the recurrence for {P 2 k) is computed as in Section 5.4, which 
gives 

{l-q)q'^z \ q3fc-i('j_g)2^2 


P2k+2 = ( 1 ~ 


P2k 


fc-ll ^'^k-2- 


(1 -I- q'=)(l + J (1-1- g'=)2(l — g2fc-i-i)(x — q'^k-iy 

The sequence Hk is defined as the numerator of the evaluation of (10), namely 

Pk{qz)Qk{z) - Pk{z)Qk{qz) „ ^ ^ \ n t \n f \ 

----- Pk(z)Qk[qz) - Qk{z)Qk[qz). 

{q-l)z 

Next, we compute a linear dependency between P[ 2 k, H 2 k+ 2 , ■ ■ ■, which is still of or¬ 
der 4 (but significantly bigger than its limit as g —5> 1). The g-analogue of reduction 
of order then proves 

g3'=+2(l — qYz'^ 


H2k+2 = —- 


TH2k. 


(1 -I- g'=+l)2(l — g.2fe-|-1^2 

which concludes the proof, providing with a generalization of the expression for exp 
that is recovered by letting g 1. 

Using the same steps leads to an automatic proof of Heine’s g-analogue of Gauss’s 
continued fraction [10, 19.2.1]: the g-hypergeometric series is defined by 

{a-,q)n{b;q)n 2 ” 


291 


{a,b;c-,q;z) = ^ 


n>0 


(c;g)„ (g;g)n’ 


where (a; g)„ is the g-Pochhammer symbol 

(a; q)n = (1 - a)(l - ag) ■ • • (1 - ag""!). 
Heine’s continued fraction is 

2 (l)i{a,b;c;q;z) 


= 1 + ^ 
2 (l)iia,bq;cq;q;z) m=i 1 


a2fe-i-i — 


(l-ag'=)(cg'=-6)g'= 

(1 — cg2^)(l — cg2^+i) ’ 


fe-i 


(1 - 6g'=)(cg'= - a)q 

(1 — cg2'=“i)(l — cg2'=) 


0'2k — 
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We sketch the main steps of the computation. The g-Riccati equation is 

(1 - cfF{z)F{qz) + (1 - c){bz - c)F{qz) 

+ (1 — z)(l — az)F{z) — z(a — l)(b — c) = 0. 

The sequence Hk oi interest is therefore the numerator of the evaluation of this 
left-hand size at F{z) = Pk{z)/Qk{z)- The continued fraction being periodic of 
period 2, a recurrence for F[ 2 k (or order 4) is computed. Reduction of order yields 

H2k+2 _ (1 - ag'"+^)(l - bq’^+^)(a - cg'"+^)(6 - cq’^+^)q'^’‘+^ 
z'^H2k (1 — cg2''+i)(l — ’ 

which concludes the proof. This automates 2 more of the formulas in [10]. 

6.3. Difference Equations. The same method applies to difference equations. 
For instance, it results in one of the classical proofs [14, chap. 3] of Brouncker’s 
continued fraction for 



where T is Euler’s Gamma function. Using the functional equation r(s-|-l) = sr(s), 
it follows that h{s) satisfies 

b{s)b{s + 2) = 16/(s -I- 1)^. 

Looking for a formal power series solution in inverse powers of s (and nonnegative 
leading term) leads to a unique solution b{s) = 4/s — 2/s^ -!-••• This is then 
converted into a continued fraction expansion with coefficients (ofc) given by 

4 1 9 25 49 81 

s ’ 2s^ ’ 4s^ ’ 4s^ ’ 4s^ ’ 4s^ ’ 

from which it is easy to conjecture Uk = {2k — 3)^/(4s^) for fc > 3. The analogue 
of Flk in this context is 

Hk = (s -l- 1)^Pk{s)Pk{s -l- 2) — 16Qk{s)Qk{s + 2), 

for which the same approach as above produces a linear recurrence of order 4 which 
is not sufficient to conclude that the valuations increase. From there, reduction of 
order with Algorithm 1 yields the shorter 

_ {2k + 1)2 

^'=+1- 4s(s + 2) ’ 

which exhibits the required increase of valuations. 

The same technique has been applied to all the explicit C-fractions concerning 
the function in [10], thereby completing the experiment on this book. 

7. Conclusion 

In a simple and unified way, our approach to continued fractions recovers an 
unexpectedly large number of explicit C-fractions from the literature. One miracle 
that takes place is that in all cases, the sequence of remainder polynomials turns 
out to be hypergeometric or g-hypergeometric. We are currently exploring this 
phenomenon in more detail. 
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